library(tidyverse)
library(bayesrules)
library(bayesplot)
library(rstan)
library(rstanarm)
library(broom.mixed)
library(tidybayes)
library(forcats)
options(mc.cores = parallel::detectCores())
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ── ✔ dplyr 1.1.2 ✔ readr 2.1.4 ✔ forcats 1.0.0 ✔ stringr 1.5.0 ✔ ggplot2 3.4.2 ✔ tibble 3.2.1 ✔ lubridate 1.9.2 ✔ tidyr 1.3.0 ✔ purrr 1.0.1 ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ── ✖ dplyr::filter() masks stats::filter() ✖ dplyr::lag() masks stats::lag() ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors This is bayesplot version 1.10.0 - Online documentation and vignettes at mc-stan.org/bayesplot - bayesplot theme set to bayesplot::theme_default() * Does _not_ affect other ggplot2 plots * See ?bayesplot_theme_set for details on theme setting Loading required package: StanHeaders rstan (Version 2.21.8, GitRev: 2e1f913d3ca3) For execution on a local, multicore CPU with excess RAM we recommend calling options(mc.cores = parallel::detectCores()). To avoid recompilation of unchanged Stan programs, we recommend calling rstan_options(auto_write = TRUE) Attaching package: ‘rstan’ The following object is masked from ‘package:tidyr’: extract Loading required package: Rcpp This is rstanarm version 2.21.4 - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors! - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults. - For execution on a local, multicore CPU with excess RAM we recommend calling options(mc.cores = parallel::detectCores()) Attaching package: ‘rstanarm’ The following object is masked from ‘package:rstan’: loo
Is there a book-level? I.e. are the characteristics of the book always the same, regardless of state and challenge?
Note that there are sometimes multiple challenges per state:
book_banning %>% filter( title %in% "A Separate Reality" )
| title | book_id | author | date | year | removed | explicit | antifamily | occult | language | lgbtq | violent | state | political_value_index | median_income | hs_grad_rate | college_grad_rate |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <fct> | <chr> | <date> | <dbl> | <fct> | <fct> | <fct> | <fct> | <fct> | <fct> | <fct> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> |
| A Separate Reality | 46 | Castaneda, Carlos | 2001-05-25 | 2001 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | OR | 4 | 1274.5 | 5.538042 | 1.07627 |
| A Separate Reality | 46 | Casteneda, Carlos | 2000-06-29 | 2000 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | OR | 4 | 1274.5 | 5.538042 | 1.07627 |
book_banning %>% filter( title %in% "Always Running" )
| title | book_id | author | date | year | removed | explicit | antifamily | occult | language | lgbtq | violent | state | political_value_index | median_income | hs_grad_rate | college_grad_rate |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <fct> | <chr> | <date> | <dbl> | <fct> | <fct> | <fct> | <fct> | <fct> | <fct> | <fct> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> |
| Always Running | 118 | Rodriguez, Luis | 2004-09-01 | 2004 | 1 | 1 | 1 | 0 | 1 | 0 | 1 | CA | 7.4 | 10119 | -2.761958 | 2.57627 |
| Always Running | 118 | Rodriguez, Luis | 2004-12-05 | 2004 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | CA | 7.4 | 10119 | -2.761958 | 2.57627 |
| Always Running | 118 | Rodriguez, Luis | 2004-10-15 | 2004 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | CA | 7.4 | 10119 | -2.761958 | 2.57627 |
Are book IDs mapped to title in unique way? In other words can multiple book IDs be assigned to the same title? Can the same title be assigned to the same book ID?
There is only one book ID for a given title.
book_banning %>% count(book_id, title) %>% count( title ) %>% filter(n > 1)
| title | n |
|---|---|
| <chr> | <int> |
However there are multiple titles assigned to the same book ID:
book_banning %>% count(book_id, title) %>% count( book_id ) %>% filter(n > 1) %>% head()
| book_id | n | |
|---|---|---|
| <fct> | <int> | |
| 1 | 100 | 2 |
| 2 | 148 | 2 |
| 3 | 356 | 3 |
| 4 | 717 | 2 |
| 5 | 722 | 2 |
| 6 | 724 | 2 |
Some examples:
book_banning %>% filter( book_id==100 ) %>% select( title, book_id, author, date )
| title | book_id | author | date |
|---|---|---|---|
| <chr> | <fct> | <chr> | <date> |
| Alice on Her Way | 100 | Naylor, Phyllis Reynolds | 2007-09-08 |
| Alice On Her Way | 100 | Naylor, Phyllis Reynolds | 2009-11-09 |
book_banning %>% filter( book_id==148 ) %>% select( title, book_id, author, date )
| title | book_id | author | date |
|---|---|---|---|
| <chr> | <fct> | <chr> | <date> |
| Angus, Thongs, and Full Frontal Snogging | 148 | Rennison, Louise | 2005-04-01 |
| Angus, Thongs, and Full Frontal Snogging | 148 | Rennison, Louise | 2003-12-11 |
| Knocked Out By My Nunga-Nungas | 148 | Rennison, Louise | 2007-07-22 |
| Angus, Thongs, and Full Frontal Snogging | 148 | Rennison, Louise | 2004-04-09 |
| Angus, Thongs, and Full Frontal Snogging | 148 | Rennison, Louise | 2009-11-22 |
| Angus, Thongs, and Full Frontal Snogging | 148 | Rennison, Louise | 2006-07-15 |
book_banning %>% filter( book_id==356 ) %>% select( title, book_id, author, date )
| title | book_id | author | date |
|---|---|---|---|
| <chr> | <fct> | <chr> | <date> |
| Captain Underpants | 356 | Pilkey, Dav | 2010-08-12 |
| Captain Underpants (series) | 356 | Pilkey, Dav | 2005-06-09 |
| Captain Underpants and the Perilous Plot of Professor Poopypants | 356 | Pilkey, Dav | 2002-06-24 |
The titles are different, but seem more or less the same for the same book ID, however this would need to be investigated more closely.
Are the book characteristics explicit, antifamily, occult, language, lgbtq and violent the same for all book challenges in a given state? To check this, compute the mean of all these variables grouped by state and book id. If the mean is different from 0 or 1, then the values per book are different. This appears not to be the case, so we can speak of a book-level that is a sublevel of state-level.
book_banning %>%
select( state, book_id, explicit:violent ) %>%
group_by( state, book_id ) %>%
summarize_all( \(x) { mean(as.numeric(x)-1) } ) %>%
filter_at( vars(explicit:violent), \(x) { !(x %in% c(1,0))} )
| state | book_id | explicit | antifamily | occult | language | lgbtq | violent |
|---|---|---|---|---|---|---|---|
| <chr> | <fct> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> |
It remains to check, whether the state-level predictors are really only on the state-level. To that end, check whether the variance of any of the predictors is larger than zero:
book_banning %>%
select( state, political_value_index, median_income, hs_grad_rate, college_grad_rate ) %>%
group_by( state ) %>%
summarize_all( var ) %>%
filter_at( vars(political_value_index:college_grad_rate), \(x) {x>0})
| state | political_value_index | median_income | hs_grad_rate | college_grad_rate |
|---|---|---|---|---|
| <chr> | <dbl> | <dbl> | <dbl> | <dbl> |
Check whether peak-level predictors do not vary from peak to peak:
climbers_sub %>%
select( peak_id, height_metres, first_ascent_year ) %>%
group_by( peak_id ) %>%
summarize_all( var ) %>%
filter_at( vars(height_metres:first_ascent_year), \(x) !(x %in% c(0,NA)) )
| peak_id | height_metres | first_ascent_year |
|---|---|---|
| <fct> | <dbl> | <dbl> |
Check whether expedition-level predictors do not vary from expedition to expedition (sublevel of peak level):
climbers_sub %>%
group_by( peak_id, expedition_id ) %>%
summarize( var_count = var(count) ) %>%
filter( var_count > 0 )
`summarise()` has grouped output by 'peak_id'. You can override using the `.groups` argument.
| peak_id | expedition_id | var_count |
|---|---|---|
| <fct> | <chr> | <dbl> |
Check whether climber-level predictors do not vary from climber to climber (sublevel of expedition level):
climbers_sub %>%
select( peak_id, expedition_id, member_id, age, expedition_role ) %>%
group_by( peak_id, expedition_id, member_id ) %>%
summarize( count=n() ) %>%
filter( count>1 )
`summarise()` has grouped output by 'peak_id', 'expedition_id'. You can override using the `.groups` argument.
| peak_id | expedition_id | member_id | count |
|---|---|---|---|
| <fct> | <chr> | <fct> | <int> |
No multiple entries for the same climber, so no variance.
spotify_small <- spotify %>%
filter(artist %in% c("Beyoncé", "Camila Cabello", "Camilo",
"Frank Ocean", "J. Cole", "Kendrick Lamar")) %>%
select(artist, album_id, popularity, valence)
head( spotify_small )
| artist | album_id | popularity | valence |
|---|---|---|---|
| <fct> | <chr> | <dbl> | <dbl> |
| Beyoncé | 7dK54iZuOxXFarGhXwEXfF | 75 | 55.2 |
| Beyoncé | 39P7VD7qlg3Z0ltq60eHp7 | 77 | 47.2 |
| Beyoncé | 1gIC63gC3B7o7FfpPACZQJ | 75 | 76.0 |
| Beyoncé | 1gIC63gC3B7o7FfpPACZQJ | 75 | 76.0 |
| Beyoncé | 20E3PwDg1jaDdK9K565luD | 66 | 47.2 |
| Beyoncé | 0Zd10MKN5j9KwUST0TdBBB | 72 | 50.9 |
dim( spotify_small )
album_id is the second grouping variable. Usually, an artist produces multiple albums with multiple songs each, so we deal with a nested model. Multiple artists can participate in an album, but then it appears to be clearly defined by providing a list in the artist field.
artists
$\to$ Beyoncé, Camila Cabello, Camilo, ..
$\to$ Album 1 of Beyoncé, Album 2 of Beyoncé, .., Album 1 of Camila Cabello, ..
$\to$ Song 1 of Album 1 of Beyoncé, .., Song 10 of Album 1 of Beyoncé, .., Song 1 of Album 1 of Camila Cabello, ..
Thus a nested structure with artists $\to$ albums $\to$ songs.
options(repr.plot.width=15, repr.plot.height=10)
ggplot( spotify_small, aes(x=valence, y=popularity, color=artist) ) +
geom_point() +
geom_line( aes(group=album_id, color=artist) )
Maybe not the nicest plot: Points represent individual songs, connected points indicate songs on the same album and colors represents artists. It appears that the variability within the songs of an individual album of an artist is smaller than the overall variability of songs/albums of an artist.
where $X_{ij}$ and $Y_{ij}$ stands for the valence and popularity for song $i$ of artist $j$. $\beta_{0j}$ represents the different baseline popularities for the different artists, $\sigma_y$ the variability within the popularity of songs of an artist and $\sigma_0$ the variability in popularity between artists. $\beta_1$ represents the global slope parameter for valence.
spotify_model_1 <- stan_glmer(
popularity ~ valence + (1 | artist),
data = spotify_small, family = gaussian,
prior_intercept = normal(0, 2.5, autoscale = TRUE),
prior = normal(0, 2.5, autoscale = TRUE),
prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1),
chains = 4, iter = 5000*2, seed = 84735
)
options(repr.plot.width=15, repr.plot.height=5)
pp_check(spotify_model_1) + xlab("popularity")
The model does not represent the measured data well! Adding album-level information might reduce the inductive bias of the model.
Here, $X_{ijk}$ and $Y_{ijk}$ stand for the valence and popularity for song $i$ in album $j$ of artist $k$. $b_{0j}$ represents the different baseline popularities for the different albums of an artist, $p_{0k}$ the different artists and $\sigma_y$ the variability within the popularity of songs of a particular album of a particular artist.$\beta_1$ represents the global slope parameter for valence.
spotify_model_2 <- stan_glmer(
popularity ~ valence + (1 | album_id) + (1 | artist),
data = spotify_small, family = gaussian,
prior_intercept = normal(0, 2.5, autoscale = TRUE),
prior = normal(0, 2.5, autoscale = TRUE),
prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1),
chains = 4, iter = 5000*2, seed = 84735
)
Warning message: “There were 6 divergent transitions after warmup. See https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup to find out why this is a problem and how to eliminate them.” Warning message: “Examine the pairs() plot to diagnose sampling problems ”
MCMC diagnostics:
options(repr.plot.width=15, repr.plot.height=15)
mcmc_trace(spotify_model_2, size = 0.1)
mcmc_dens_overlay(spotify_model_2)
mcmc_acf(spotify_model_2)
neff_ratio(spotify_model_2)
rhat(spotify_model_2)
Warning message: “The `facets` argument of `facet_grid()` is deprecated as of ggplot2 2.2.0. ℹ Please use the `rows` argument instead. ℹ The deprecated feature was likely used in the bayesplot package. Please report the issue at <https://github.com/stan-dev/bayesplot/issues/>.”
Looks ok, ignore the warning, however in practice one should tweak the MCMC parameters here.
Posterior predictive check:
pp_check(spotify_model_2) + xlab("popularity")
Wow, this is much better!
Fixed effects:
tidy(spotify_model_2, effects = "fixed")
| term | estimate | std.error |
|---|---|---|
| <chr> | <dbl> | <dbl> |
| (Intercept) | 61.95850503 | 5.01287305 |
| valence | 0.06480214 | 0.03787398 |
Random effects:
tidy(spotify_model_2, effects = "ran_vals") %>%
filter( level %in% c("Kendrick_Lamar", "748dZDqSZy6aPXKcI9H80u") )
| level | group | term | estimate | std.error |
|---|---|---|---|---|
| <chr> | <chr> | <chr> | <dbl> | <dbl> |
| 748dZDqSZy6aPXKcI9H80u | album_id | (Intercept) | 19.02559 | 7.044108 |
| Kendrick_Lamar | artist | (Intercept) | -16.42305 | 7.635712 |
spotify_small dataset: popularity = 62.0 + 0.065 * valenceKendrick Lamar has in general a worse average popularity than the average of our sample. However his Album “good kid, m.A.A.d city” was so popular that a song in it has a larger popularity than the average of our sample.
tidy(spotify_model_2, effects = "ran_vals") %>%
filter( group == "artist" ) %>%
arrange( desc(estimate) )
| level | group | term | estimate | std.error |
|---|---|---|---|---|
| <chr> | <chr> | <chr> | <dbl> | <dbl> |
| Camilo | artist | (Intercept) | 6.298670 | 7.279473 |
| Camila_Cabello | artist | (Intercept) | 3.748713 | 5.899080 |
| J._Cole | artist | (Intercept) | 3.601396 | 6.845808 |
| Beyoncé | artist | (Intercept) | 2.820226 | 5.765683 |
| Frank_Ocean | artist | (Intercept) | 1.791417 | 5.608984 |
| Kendrick_Lamar | artist | (Intercept) | -16.423055 | 7.635712 |
Camilo.
Most popular album:
tidy(spotify_model_2, effects = "ran_vals") %>%
filter( group == "album_id" ) %>%
arrange( desc(estimate) ) %>%
head()
| level | group | term | estimate | std.error |
|---|---|---|---|---|
| <chr> | <chr> | <chr> | <dbl> | <dbl> |
| 4eLPsYPBmXABThSJ821sqY | album_id | (Intercept) | 30.25010 | 6.429640 |
| 3pLdWdkj83EYfDN6H2N8MR | album_id | (Intercept) | 26.84868 | 7.875206 |
| 748dZDqSZy6aPXKcI9H80u | album_id | (Intercept) | 19.02559 | 7.044108 |
| 6PBZN8cbwkqm1ERj2BGXJ1 | album_id | (Intercept) | 18.18493 | 7.323289 |
| 3XzSOIE6zGLliuqsVGLmUc | album_id | (Intercept) | 16.67567 | 7.243112 |
| 3Vsbl0diFGw8HNSjG8ue9m | album_id | (Intercept) | 16.50975 | 5.246488 |
spotify %>%
filter( album_id == "4eLPsYPBmXABThSJ821sqY" ) %>%
select( artist, album_name ) %>%
first()
| artist | album_name |
|---|---|
| <fct> | <chr> |
| Kendrick Lamar | DAMN. |
tidy(spotify_model_2, effects = "ran_pars") %>%
mutate( squared = estimate^2 ) %>%
mutate( proportion = squared/sum(squared) ) %>%
select( -squared )
| term | group | estimate | proportion |
|---|---|---|---|
| <chr> | <chr> | <dbl> | <dbl> |
| sd_(Intercept).album_id | album_id | 17.299246 | 0.65699327 |
| sd_(Intercept).artist | artist | 11.029960 | 0.26708805 |
| sd_Observation.Residual | Residual | 5.880592 | 0.07591869 |
$\sigma_y = 5.9, \sigma_b = 17.3, \sigma_p = 11.0$ ($\sigma_b$ represents variability between different albums of an artist and $\sigma_p$ variability between artists).
The largest proportion of variability in popularity (65%) comes from $\sigma_b$, i.e. the albums of an artist, followed by $\sigma_p$, the variability of popularity between artists. The variability of the popularity of songs within an album is interestingly the smallest. It appears thus that most artists do not consistently produce good albums, but rather a diverse range of popular and unpopular albums. The variance in the individual songs of an album is rather small (7.6%) compared with the other effects, i.e. if one song in an album is popular, usually the rest is also popular and the other way round.
bwc <- big_word_club %>%
mutate(treat = as.factor(treat))
Interestingly, treat is a school-level predictor and the same for all students at a particular school.
bwc %>%
group_by( school_id ) %>%
summarize(
var_private_school = var(private_school),
var_lunch = var(free_reduced_lunch),
var_treat = var(as.numeric(treat)-1),
var_age_months = var(age_months),
var_esl_observed = var(esl_observed)
) %>%
mutate_if( is.factor, \(x) as.numeric(x)-1 ) %>%
summarize_all( max ) %>%
select(-school_id)
| var_private_school | var_lunch | var_treat | var_age_months | var_esl_observed |
|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> |
| 0 | 0 | 0 | 174.8603 | 0.2614379 |
No variance among the school-level predictors, variance among the student-level predictors.
How many tests per student?
bwc %>%
count( participant_id ) %>%
filter( n>1 )
| participant_id | n |
|---|---|
| <int> | <int> |
Only one! Not a nested grouping structure.
Use a normal hierarchical regression model with random intercepts. For simplicity, do not use random slopes.
Here, $Y_{ij}$ stands for percent change in score for student $i$ at school $j$. $b_{0j}$ represents the different baseline percent changes for the different schools and $\sigma_y$ the variability within the different students at an individual school. $\beta_1$, $\beta_2$ and $\beta_3$ are the global slope parameters for the school-level predictors treat, private_school and free_reduced_lunch (denoted as $U_{j1}$, $U_{j2}$ and $U_{j3}$). $\beta_4$ and $\beta_5$ are the slopes for the student-level predictors age_months and esl_observed (denoted as $X_{ij1}$ and $X_{ij2}$).
bwc_model <- stan_glmer(
score_pct_change ~ age_months + esl_observed + treat + private_school +
free_reduced_lunch + (1 | school_id),
data = bwc, family = gaussian,
prior_intercept = normal(0, 2.5, autoscale = TRUE),
prior = normal(0, 2.5, autoscale = TRUE),
prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1),
chains = 4, iter = 5000*2, seed = 84735,
)
Diagnostics:
options(repr.plot.width=15, repr.plot.height=15)
mcmc_trace(bwc_model, size = 0.1)
mcmc_dens_overlay(bwc_model)
mcmc_acf(bwc_model)
neff_ratio(bwc_model)
rhat(bwc_model)
Looks good!
Posterior predictive check for entire model:
options(repr.plot.width=15, repr.plot.height=5)
pp_check(bwc_model) + xlab("score percent change")
Reasonable, although the normal does not capture the data distribution perfectly.
tidy(bwc_model, effects = "fixed", conf.int = TRUE, conf.level = 0.80)
| term | estimate | std.error | conf.low | conf.high |
|---|---|---|---|---|
| <chr> | <dbl> | <dbl> | <dbl> | <dbl> |
| (Intercept) | 17.61253771 | 6.46763770 | 9.43693441 | 25.81506041 |
| age_months | -0.17189181 | 0.08883075 | -0.28411372 | -0.05867890 |
| esl_observed | 0.73754942 | 2.80741039 | -2.83367186 | 4.34695532 |
| treat1 | -0.88277744 | 1.63783580 | -2.97946784 | 1.23215486 |
| private_school | 0.89210981 | 3.07991282 | -3.15180117 | 4.83356378 |
| free_reduced_lunch | 0.01817009 | 0.03533995 | -0.02757163 | 0.06304653 |
Only age_months appears to be significant at the 80% level, with a large uncertainty. It appears that with each month of age, score_pct_change is reduced by 0.17.
tidy(bwc_model, effects = "ran_pars") %>%
mutate( squared=estimate^2, proportion=squared/sum(squared) ) %>%
select( -squared )
| term | group | estimate | proportion |
|---|---|---|---|
| <chr> | <chr> | <dbl> | <dbl> |
| sd_(Intercept).school_id | school_id | 2.509992 | 0.02033181 |
| sd_Observation.Residual | Residual | 17.423032 | 0.97966819 |
$\sigma_y=17.4$, $\sigma_b=2.5$. The variability between the students of a school clearly dominates the overall variability.
school_level_effects <- tidy(bwc_model, effects = "ran_vals", conf.int = TRUE, conf.level = 0.80)
head( school_level_effects )
| level | group | term | estimate | std.error | conf.low | conf.high |
|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> |
| 1 | school_id | (Intercept) | -0.2128394 | 1.595665 | -3.2222516 | 1.9455620 |
| 2 | school_id | (Intercept) | 0.9269347 | 1.816862 | -0.8907846 | 4.6999950 |
| 3 | school_id | (Intercept) | 0.8832702 | 1.820799 | -0.9400862 | 4.6748250 |
| 4 | school_id | (Intercept) | 0.1742617 | 1.596676 | -2.0269883 | 3.0610136 |
| 5 | school_id | (Intercept) | 0.3163362 | 1.607116 | -1.6833285 | 3.4417448 |
| 6 | school_id | (Intercept) | -1.2107841 | 2.042551 | -5.5164243 | 0.6814634 |
school_level_effects %>%
filter( level == "30" )
| level | group | term | estimate | std.error | conf.low | conf.high |
|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> |
| 30 | school_id | (Intercept) | -1.591776 | 2.353159 | -6.042635 | 0.3968267 |
Posterior median model: score_pct_change = 17.6 - 1.6 - 0.17 age_months = 16.0 - 0.17 age_months (reporting only significant coefficients)
school_level_effects %>%
filter( level == "47" )
| level | group | term | estimate | std.error | conf.low | conf.high |
|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> |
| 47 | school_id | (Intercept) | 1.354891 | 2.144308 | -0.5696326 | 5.652695 |
Posterior median model: score_pct_change = 17.6 + 1.4 - 0.17 age_months = 19.0 - 0.17 age_months (reporting only significant coefficients)
Posterior median model: score_pct_change = 17.6 - 0.17 age_months + 0.018 * 95 = 19.3 - 0.17 age_months
(using the non-significant contribution of free_reduced_lunch)
Posterior median model: score_pct_change = 17.6 - 0.17 age_months + 0.018 * 10 = 17.8 - 0.17 age_months
(using the non-significant contribution of free_reduced_lunch)
None at a significant level. Maybe with more data private_school, maybe also the availability of free or reduced lunch, however I do not see this directly.
The younger the student, the greater the vocabular improvement (at a significant level, yet with quite some uncertainty). Maybe with more data it also helps being an ESL student (to have more change / maybe a steeper learning curve).